nz=240
nx=120
Lz=1.4
Lx=5.
pi=3.141592

print, 'setting cutoff wavenumbers in z direction ...'
kzmin=2*pi/Lz
kzmax=nz*kzmin

kxmin=2*pi/Lx
kxmax=nx*kxmin

kx=[indgen(nx/2+1),-reverse(indgen(nx/2-1)+1)]*kxmin
kz=[indgen(nz/2+1),-reverse(indgen(nz/2-1)+1)]*kzmin
ks=indgen(nz/2)

len=size(u)
nt=len[4]

umod = fltarr(nx,nz,nt)
wmod = fltarr(nx,nz,nt)

for it=0,nt-1 do begin
 for ix=0,nx-1 do begin
  for iz=0,nz-1 do begin
   umod[ix,iz,it] = u[ix,0,iz,it]
   wmod[ix,iz,it] = w[ix,0,iz,it] 
  endfor
 endfor
endfor

ftu = fltarr(nx,nz,nt)
ftw = fltarr(nx,nz,nt)
E = fltarr(n_elements(ks))
;Eprom = fltarr(nx,nz)

for it=0,nt-1 do begin
 print, 'computing 2D fft for temporal slice', it
 ftu[*,*,it] = fft(umod[*,*,it],/double)
 ftw[*,*,it] = fft(wmod[*,*,it],/double)
endfor

init=10
print, 'computing energy amplitudes ...'
for it=nt-init,nt-1 do begin
 for ix=0,nx-1 do begin
  for iz=0,nz-1 do begin
   k=round(sqrt(kx[ix]^2+kz[iz]^2))
   if(k lt n_elements(ks)) then begin
    E[k]=E[k]+abs(ftu[ix,iz,it])^2+abs(ftw[ix,iz,it])^2
   endif
  endfor
 endfor
endfor

Er=E/double(nx*init)

kzlist=fltarr(nz)
for iz=0, nz-1 do begin
 kzlist[iz]=kzmin*(iz+1)
endfor

print,  'preparing plot of energy spectra....'

!x.margin=[9,3]
!p.charsize=1.5

names=['k!U-5/3!N','k!U-3!N']
;names=['k^-(5/3)','k^-3']
xpos=[90,30]
ypos=[5e-7,5e-6]

intercept1=6e-4
intercept2=9e-2
yrmax=max(Er)
yrmin=9e-9
plot, Er, /xlog, /ylog,xr=[2,150], xtitle='Wavenumber k',ytitle='E(k)',xstyle=1,yr=[yrmin,yrmax]
oplot, intercept1*ks^(-5./3.),linestyle=2, color=160
oplot, intercept2*ks^(-3.), linestyle=2, color='FF0000'x
for ii=0,1 do xyouts, xpos[ii], ypos[ii], names[ii]

end
